We analysed 1434 flights of 50 individuals of blue tit, great tit, and marsh tit as baseline data We found significant difference in initial velocity and in-flight acceleration between species . Great tits showed higher initial velocity and lower in-flight acceleration than blue tits. Marsh tits did not differ in flight performance with other species. Weather conditions did not show significant effect to flight performance.
library(dplyr)
library(stringr)
library(readr)
library(fuzzyjoin)
library(tidyverse)
library(ggplot2)
library(lubridate)
library(ggpubr)
library(rstatix)
library(coin)
library(scales)
library(gridExtra)
library(ggthemes)
library(ggbreak)
library(lme4)
library(reshape2)
library(olsrr)
library(gvlma)
library(PerformanceAnalytics)
library(usethis)
library(devtools)
library(redres)
library(EnvStats)
library(lmtest)
library(lmerTest)
library(AICcmodavg)
library(emmeans)
library(MuMIn)
library(sjPlot)
#in dept laptop
#setwd("C:\\Users\\jesu4036\\OneDrive - Nexus365\\Documents\\R\\WythamWoods")
#load all dates data
#all_dates3 <- read.delim(file = ".\\all_dates3.txt", sep = ",")
#in the kyubird-laptop
all_dates3 <- read.delim(file = "C:\\Users\\kmh\\Desktop\\DATA\\all_dates3.txt", sep = ",") %>% distinct()
all_dates3$datetime.y <- as.POSIXct(all_dates3$datetime.y, format = "%Y-%m-%d %H:%M:%S")
colnames(all_dates3)
## [1] "Ring" "RFID" "datetime.y" "temperature" "firstbreak1"
## [6] "firstbreak2" "firstbreak3" "lastbreak1" "lastbreak2" "lastbreak3"
## [11] "t1" "t2" "v1" "v2" "avg.v"
## [16] "acc" "section2" "section3" "section4" "length1"
## [21] "length2" "tunnel" "species" "age" "Spec"
## [26] "Sex" "Date" "Place" "Site" "Wing"
## [31] "Wt" "Init" "Tarsus"
all_dates3 <- all_dates3 %>%
select(-Spec)
#load climate data
#in dept laptop
#climate <- read.delim(file = ".\\climate_filled.csv", sep = ',')
#in the kyubird laptop
climate <- read.delim(file = "C:\\Users\\kmh\\Desktop\\DATA\\climate_filled.csv", sep = ',')
climate$datetime <- as.POSIXct(climate$datetime, format = "%Y-%m-%d %H:%M:%S")
#in the dept laptop
#caught.hour.day <- read.delim(file = ".\\caught.hour.day.csv", sep = ",")
#in the kyubird laptop
caught.hour.day <- read.delim(file = "C:\\Users\\kmh\\Desktop\\DATA\\caught.hour.day.csv", sep = ',')
all_dates3$day <- day(all_dates3$datetime.y)
all_dates3$hour <- hour(all_dates3$datetime.y)
all_dates3$roundate <- round_date(all_dates3$datetime.y, "10 mins")
all_dates3.climate <- merge(x = all_dates3, y = climate, by.x = "roundate", by.y = "datetime", all.x = TRUE)
flight.data <- all_dates3.climate %>%
filter(day %in% c(15:24))
nat.var <- all_dates3.climate %>%
filter(day %in% c(15:19))
capt.eff <- all_dates3.climate %>%
filter(day %in% c(20:24))
write.csv(nat.var, “C:\Users\kmh\Desktop\DATA\nat.var.csv”)
dim(flight.data) #total number of flight
## [1] 3451 38
table(flight.data$day)
##
## 15 16 17 18 19 20 21 22 23 24
## 241 231 54 528 380 485 616 318 504 94
no.of.individuals <- rowSums(table(flight.data$species, flight.data$Ring)!=0)
no.of.flights <- table(flight.data$species)
raw_sum <- rbind(no.of.individuals, no.of.flights)
raw_sum
## bluti greti marti
## no.of.individuals 40 16 5
## no.of.flights 2617 513 321
dim(nat.var) #total number of flight
## [1] 1434 38
table(nat.var$day)
##
## 15 16 17 18 19
## 241 231 54 528 380
no.of.individuals <- rowSums(table(nat.var$species, nat.var$Ring)!=0)
no.of.flights <- table(nat.var$species)
nat.var_sum <- rbind(no.of.individuals, no.of.flights)
nat.var_sum
## bluti greti marti
## no.of.individuals 33 12 4
## no.of.flights 1172 148 114
breeding <- cbind(nat.var$Ring, nat.var$RFID, nat.var$species, nat.var$age, nat.var$Sex) %>% as.data.frame() %>%
distinct()
colnames(breeding) <- c("Ring", "RFID", "species", "age", "Sex")
write.csv(breeding, “\breeding.csv”,row.names = FALSE )
#find the ring number of those without Sex in Great Tit & Blue Tit
for( i in 1:nrow(nat.var)){
if(nat.var$Ring[i] == "ATR8387"){
nat.var$Sex[i] <- "F"
}
else if(nat.var$Ring[i] == "ATR8383"){
nat.var$Sex[i] <- "F"
}
else if(nat.var$Ring[i] == "PF72506"){
nat.var$Sex[i] <- "M"
}
else{}
}
#Sex unknown for ATR8380 and ATR5368 both in morphometrics and brood data
nosex <- nat.var %>%
filter(species != "marti") %>%
filter(Sex == "") %>%
select(Ring) %>%
distinct()
numbers <- nat.var %>%
count(species, age, Sex, Ring)
individual.age <- numbers %>%
count(species, age)
individual.sex <- numbers %>%
count(species, Sex)
numbers.age <- numbers %>%
group_by(species, age) %>%
summarise(mean = mean(n), sd = sd(n)) %>%
na.omit()
## `summarise()` has grouped output by 'species'. You can override using the
## `.groups` argument.
numbers.sex <- numbers %>%
group_by(species, Sex) %>%
summarise(mean = mean(n), sd = sd(n)) %>%
na.omit()
## `summarise()` has grouped output by 'species'. You can override using the
## `.groups` argument.
numbers.age <- numbers.age %>%
add_column(individual.age$n)
numbers.sex <- numbers.sex %>%
add_column(individual.sex$n)
numbers.tot <- numbers %>%
group_by(species) %>%
summarise(mean = mean(n), sd = sd(n)) %>%
na.omit()
numbers.age
## # A tibble: 5 × 5
## # Groups: species [3]
## species age mean sd `individual.age$n`
## <chr> <int> <dbl> <dbl> <int>
## 1 bluti 5 33.6 24.6 22
## 2 bluti 6 36 20.0 12
## 3 greti 5 11.4 11.5 10
## 4 greti 6 17 22.6 2
## 5 marti 5 28.5 11.0 4
numbers.sex
## # A tibble: 6 × 5
## # Groups: species [3]
## species Sex mean sd `individual.sex$n`
## <chr> <chr> <dbl> <dbl> <int>
## 1 bluti "" 54.5 17.7 2
## 2 bluti "F" 33.7 18.9 15
## 3 bluti "M" 32.8 26.2 17
## 4 greti "F" 15.8 17.4 4
## 5 greti "M" 10.6 10.5 8
## 6 marti "" 28.5 11.0 4
numbers.tot
## # A tibble: 3 × 3
## species mean sd
## <chr> <dbl> <dbl>
## 1 bluti 34.5 22.8
## 2 greti 12.3 12.6
## 3 marti 28.5 11.0
bluti <- numbers %>%
filter(species == "bluti")
greti <- numbers %>%
filter(species == "greti")
marti <- numbers %>%
filter(species == "marti")
bluti.F <- numbers %>%
filter(species == "bluti" & Sex == "F")
bluti.M <- numbers %>%
filter(species == "bluti" & Sex == "M")
bluti.5 <- numbers %>%
filter(species == "bluti" & age == 5)
bluti.6 <- numbers %>%
filter(species == "bluti" & age == 6)
greti.F <- numbers %>%
filter(species == "greti" & Sex == "F")
greti.M <- numbers %>%
filter(species == "greti" & Sex == "M")
greti.5 <- numbers %>%
filter(species == "greti" & age == 5)
greti.6 <- numbers %>%
filter(species == "greti" & age == 6)
t.test(bluti.F$n, bluti.M$n, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
##
## Welch Two Sample t-test
##
## data: bluti.F$n and bluti.M$n
## t = 0.10527, df = 28.948, p-value = 0.9169
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -15.53856 17.22483
## sample estimates:
## mean of x mean of y
## 33.66667 32.82353
t.test(bluti.5$n, bluti.6$n, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
##
## Welch Two Sample t-test
##
## data: bluti.5$n and bluti.6$n
## t = -0.30328, df = 27.045, p-value = 0.764
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -18.35356 13.62629
## sample estimates:
## mean of x mean of y
## 33.63636 36.00000
t.test(greti.F$n, greti.M$n, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
##
## Welch Two Sample t-test
##
## data: greti.F$n and greti.M$n
## t = 0.5423, df = 4.1323, p-value = 0.6155
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -20.7855 31.0355
## sample estimates:
## mean of x mean of y
## 15.750 10.625
t.test(greti.5$n, greti.6$n, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
##
## Welch Two Sample t-test
##
## data: greti.5$n and greti.6$n
## t = -0.34135, df = 1.1049, p-value = 0.7863
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -172.6219 161.4219
## sample estimates:
## mean of x mean of y
## 11.4 17.0
t.test(bluti$n, greti$n, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
##
## Welch Two Sample t-test
##
## data: bluti$n and greti$n
## t = 4.146, df = 35.353, p-value = 0.0002014
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 11.30160 32.97291
## sample estimates:
## mean of x mean of y
## 34.47059 12.33333
t.test(bluti$n, marti$n, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
##
## Welch Two Sample t-test
##
## data: bluti$n and marti$n
## t = 0.88326, df = 6.6185, p-value = 0.408
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -10.20223 22.14341
## sample estimates:
## mean of x mean of y
## 34.47059 28.50000
t.test(greti$n, marti$n, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
##
## Welch Two Sample t-test
##
## data: greti$n and marti$n
## t = -2.447, df = 5.8748, p-value = 0.05085
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -32.41668686 0.08335352
## sample estimates:
## mean of x mean of y
## 12.33333 28.50000
#overall climate data
#find how many weather values are missing
#missing data that do not make to the final dataset.
View(flight.data[!complete.cases(flight.data),])
#mean weather
mean(flight.data$air.temp, na.rm = TRUE)
## [1] 9.096442
mean(flight.data$headwind, na.rm = TRUE) #relative to the birds flying out of the tunnel
## [1] -1.428632
mean(flight.data$humidity, na.rm = TRUE)
## [1] 80.46904
#dataset for baseline variation
climate.nat.var <- nat.var %>%
filter(day %in% c(15:19)) %>%
select(headwind, air.temp, humidity)
#not using caught.hour.day which is a reduced dataset to just answer the question if the climate conditions were different at nat.var vs. capt.eff
#otherwise everything is significant probably the sample size is so vastly different
climate.capt.eff <- caught.hour.day %>%
select(headwind, air.temp, humidity)
#mean of nat.var. # 2 for columns
apply(climate.nat.var, 2,mean)
## headwind air.temp humidity
## -3.169797 10.976903 78.366511
#sd of nat.var
apply(climate.nat.var, 2, sd)
## headwind air.temp humidity
## 1.040167 1.553028 9.631564
#mean of capt.eff
apply(climate.capt.eff, 2, mean, na.rm = TRUE)
## headwind air.temp humidity
## 2.300897 6.819277 72.716639
#mean of capt.eff
apply(climate.capt.eff, 2, sd, na.rm = TRUE)
## headwind air.temp humidity
## 2.1071941 0.9258642 5.7299068
climate.all <- nat.var %>%
filter(day %in% c(15:24)) %>%
select(headwind, air.temp, humidity)
apply(climate.all, 2, mean)
## headwind air.temp humidity
## -3.169797 10.976903 78.366511
apply(climate.all, 2, sd)
## headwind air.temp humidity
## 1.040167 1.553028 9.631564
caught.byday <- caught.hour.day %>%
group_by(cgt_day) %>%
summarise(mean.hdw = mean(headwind, na.rm = FALSE),
mean.art = mean(air.temp, na.rm = FALSE),
mean.hmd = mean(humidity, na.rm = FALSE),
sd.hdw = sd(headwind, na.rm = FALSE),
sd.art = sd(air.temp, na.rm = FALSE),
sd.hmd = sd(humidity, na.rm = FALSE),
max.hdw = max(headwind, na.rm = FALSE),
min.hdw = min(headwind, na.rm = FALSE)
) %>%
na.omit()
caught.byday
## # A tibble: 2 × 9
## cgt_day mean.hdw mean.art mean.hmd sd.hdw sd.art sd.hmd max.hdw min.hdw
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 0.156 7.58 75.8 1.53 0.774 5.46 4.83 -1.05
## 2 1 3.71 6.32 70.7 0.894 0.633 4.98 5.47 1.99
day0 <- caught.hour.day %>%
filter(cgt_day == 0)
day1 <- caught.hour.day %>%
filter(cgt_day == 1)
t.test(day0$headwind, day1$headwind, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
##
## Welch Two Sample t-test
##
## data: day0$headwind and day1$headwind
## t = -18.211, df = 107.5, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.942435 -3.168404
## sample estimates:
## mean of x mean of y
## 0.1563588 3.7117780
t.test(day0$air.temp, day1$air.temp, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
##
## Welch Two Sample t-test
##
## data: day0$air.temp and day1$air.temp
## t = 11.723, df = 136.2, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 1.045330 1.469544
## sample estimates:
## mean of x mean of y
## 7.577731 6.320294
t.test(day0$humidity, day1$humidity, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
##
## Welch Two Sample t-test
##
## data: day0$humidity and day1$humidity
## t = 6.4736, df = 148.01, p-value = 1.323e-09
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 3.527083 6.626559
## sample estimates:
## mean of x mean of y
## 75.77885 70.70203
t.test(climate.nat.var$headwind, climate.capt.eff$headwind, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
##
## Welch Two Sample t-test
##
## data: climate.nat.var$headwind and climate.capt.eff$headwind
## t = -35.132, df = 200.24, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -5.777751 -5.163637
## sample estimates:
## mean of x mean of y
## -3.169797 2.300897
t.test(climate.nat.var$air.temp, climate.capt.eff$air.temp, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
##
## Welch Two Sample t-test
##
## data: climate.nat.var$air.temp and climate.capt.eff$air.temp
## t = 52.727, df = 347.02, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 4.002540 4.312713
## sample estimates:
## mean of x mean of y
## 10.976903 6.819277
t.test(climate.nat.var$humidity, climate.capt.eff$humidity, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
##
## Welch Two Sample t-test
##
## data: climate.nat.var$humidity and climate.capt.eff$humidity
## t = 11.571, df = 347.77, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 4.689545 6.610200
## sample estimates:
## mean of x mean of y
## 78.36651 72.71664
raw.flight.metrics <- nat.var %>%
group_by(species) %>%
summarise(mean.v1 = mean(v1, na.rm = FALSE),
mean.v2 = mean(v2, na.rm = FALSE),
mean.acc = mean(acc, na.rm = FALSE),
mean.avg.v = mean(avg.v, na.rm = FALSE),
mean.t1 = mean(t1, na.rm = FALSE),
mean.t2 = mean(t2, na.rm = FALSE),
mean.t3 = mean((t1+t2), na.rm = FALSE),
sd.v1 = sd(v1, na.rm = FALSE),
sd.v2 = sd(v2, na.rm = FALSE),
sd.acc = sd(acc, na.rm = FALSE),
sd.avg.v = sd(avg.v, na.rm = FALSE),
sd.t1 = sd(t1, na.rm = FALSE),
sd.t2 = sd(t2, na.rm = FALSE),
sd.t3 = sd((t1+t2), na.rm = FALSE)
) %>%
na.omit()
raw.flight.metrics
## # A tibble: 3 × 15
## species mean.v1 mean.v2 mean.acc mean.avg.v mean.t1 mean.t2 mean.t3 sd.v1
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 bluti 2.53 2.94 3.89 2.71 107. 105. 212. 0.451
## 2 greti 2.84 2.87 -0.0160 2.83 97.4 108. 206. 0.624
## 3 marti 2.62 2.85 1.96 2.71 106. 108. 214. 0.511
## # ℹ 6 more variables: sd.v2 <dbl>, sd.acc <dbl>, sd.avg.v <dbl>, sd.t1 <dbl>,
## # sd.t2 <dbl>, sd.t3 <dbl>
bluti.raw <- nat.var %>%
filter(species == "bluti") %>%
mutate(t3 = t1 + t2)
greti.raw <- nat.var %>%
filter(species == "greti") %>%
mutate(t3 = t1 + t2)
t.test(bluti.raw$t1, greti.raw$t1, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
##
## Welch Two Sample t-test
##
## data: bluti.raw$t1 and greti.raw$t1
## t = 3.943, df = 192.45, p-value = 0.0001126
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 4.857541 14.581072
## sample estimates:
## mean of x mean of y
## 107.13823 97.41892
t.test(bluti.raw$t2, greti.raw$t2, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
##
## Welch Two Sample t-test
##
## data: bluti.raw$t2 and greti.raw$t2
## t = -1.6142, df = 207.97, p-value = 0.108
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -7.0957960 0.7070403
## sample estimates:
## mean of x mean of y
## 104.9881 108.1824
t.test(bluti.raw$t3, greti.raw$t3, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
##
## Welch Two Sample t-test
##
## data: bluti.raw$t3 and greti.raw$t3
## t = 1.6545, df = 188.86, p-value = 0.09969
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.25462 14.30448
## sample estimates:
## mean of x mean of y
## 212.1263 205.6014
write.csv(raw.flight.metrics, “\raw.flight.metrics.csv”,row.names = FALSE )
day effect not included in linear or linear mixed models because
there was no difference among days until mist-netting happened.
We chose linear mixed model to account for individual variation in flight speed. We chose headwind from climate factors by its relationship with flight and correlation with air temperature and humidity.
## [1] "roundate" "Ring" "RFID" "datetime.y" "temperature"
## [6] "firstbreak1" "firstbreak2" "firstbreak3" "lastbreak1" "lastbreak2"
## [11] "lastbreak3" "t1" "t2" "v1" "v2"
## [16] "avg.v" "acc" "section2" "section3" "section4"
## [21] "length1" "length2" "tunnel" "species" "age"
## [26] "Sex" "Date" "Place" "Site" "Wing"
## [31] "Wt" "Init" "Tarsus" "day" "hour"
## [36] "air.temp" "headwind" "humidity"
lmer.v1 = lmer(v1 ~ species + Sex + age + headwind + (1|Ring), data = nat.var)
lmer.v2 = lmer(v2 ~species + Sex + age + headwind + (1|Ring), data = nat.var)
lmer.avg.v = lmer(avg.v ~ species + Sex + age + headwind + (1|Ring), data = nat.var)
lmer.acc = lmer(acc ~ species + Sex + age + headwind + (1|Ring), data = nat.var)
lmer.v1 = lmer(v1 ~ species + Sex + age + air.temp + (1|Ring), data = nat.var)
lmer.v2 = lmer(v2 ~species + Sex + age + air.temp + (1|Ring), data = nat.var)
lmer.avg.v = lmer(avg.v ~ species + Sex + age + air.temp + (1|Ring), data = nat.var)
lmer.acc = lmer(acc ~ species + Sex + age + air.temp + (1|Ring), data = nat.var)
We then identified 15 outliers using Rosner’s Test(1983) and prepared a dataset excluding outliers to see the sensitivity of the models to the outliers. All the outliers were included in the dataset in the end as we could not confirm that these outliers were from measurement errors.
test.v1 <- rosnerTest(nat.var$v1, k=3)
test.v1$all.stats
test.v2 <- rosnerTest(nat.var$v2, k=10)
test.v2$all.stats
test.avg.v<- rosnerTest(nat.var$avg.v, k=13)
test.avg.v$all.stats
test.acc <- rosnerTest(nat.var$acc, k=3)
test.acc$all.stats
#collecting outlier row numbers
out.v1 <- test.v1$all.stats %>%
filter(Outlier == "TRUE") %>%
select(Obs.Num)
out.v2 <- test.v2$all.stats %>%
filter(Outlier == "TRUE") %>%
select(Obs.Num)
out.avg.v <- test.avg.v$all.stats %>%
filter(Outlier == "TRUE") %>%
select(Obs.Num)
out.acc <- test.acc$all.stats %>%
filter(Outlier == "TRUE") %>%
select(Obs.Num)
outliers <- rbind(out.v1, out.v2, out.avg.v, out.acc) %>% distinct()
outliers #15 outliers
## Obs.Num
## 1 1014
## 2 243
## 3 1318
## 4 951
## 5 1088
## 6 1177
## 7 27
## 8 1333
## 9 351
## 10 861
## 11 1222
## 12 1257
## 13 1207
## 14 1315
## 15 195
nat.var.out <- nat.var[-c(outliers$Obs.Num),]
dim(nat.var.out) #1419 observations
## [1] 1419 38
# Computes the default residuals (raw conditional)
raw_cond.v1 <- compute_redres(lmer.v1)
## Loading required namespace: testthat
raw_cond.v2 <- compute_redres(lmer.v2)
raw_cond.avg.v <- compute_redres(lmer.avg.v)
raw_cond.acc <- compute_redres(lmer.acc)
shapiro.test(raw_cond.v1) #p <0.05
##
## Shapiro-Wilk normality test
##
## data: raw_cond.v1
## W = 0.9694, p-value < 2.2e-16
shapiro.test(raw_cond.v2) #p <0.05
##
## Shapiro-Wilk normality test
##
## data: raw_cond.v2
## W = 0.89926, p-value < 2.2e-16
shapiro.test(raw_cond.avg.v) #p <0.05
##
## Shapiro-Wilk normality test
##
## data: raw_cond.avg.v
## W = 0.89764, p-value < 2.2e-16
shapiro.test(raw_cond.acc) #p <0.05
##
## Shapiro-Wilk normality test
##
## data: raw_cond.acc
## W = 0.99357, p-value = 7.024e-06
plot(nat.var$headwind, nat.var$v1) #linearity
plot(nat.var$air.temp, nat.var$v1) #linearity
hist(raw_cond.v1) #normality in histogram
plot_redres(lmer.v1) #heteroskedasticity
plot_resqq(lmer.v1) #normality of residuals
plot_ranef(lmer.v1) #normality of random effects residuals
plot(nat.var$headwind, nat.var$v2)
plot(nat.var$air.temp, nat.var$v2) #wider range
hist(raw_cond.v2) #normality in histogram
plot_redres(lmer.v2) #heteroskedasticity
plot_resqq(lmer.v2) #normality of residuals. not great..
plot_ranef(lmer.v2) #normality of random effects residuals
plot(nat.var$headwind, nat.var$avg.v)
plot(nat.var$air.temp, nat.var$avg.v)
hist(raw_cond.avg.v) #normality in histogram
plot_redres(lmer.avg.v) #heteroskedasticity
plot_resqq(lmer.avg.v) #normality of residuals. #not great..
plot_ranef(lmer.avg.v) #normality of random effects residuals
plot(nat.var$headwind, nat.var$acc)
plot(nat.var$air.temp, nat.var$acc) #better linearity
hist(raw_cond.acc) #normality in histogram
plot_redres(lmer.acc) #heteroskedasticity
plot_resqq(lmer.acc) #normality of residuals. great.
plot_ranef(lmer.acc) #normality of random effects residuals
library(lmerTest)
Anova(lmer.v1, type = 3, test.statistic = "F") #Kenward-Rodger df
## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
##
## Response: v1
## F Df Df.res Pr(>F)
## (Intercept) 29.7686 1 40.38 2.671e-06 ***
## species 4.7090 2 39.58 0.01464 *
## Sex 0.9514 2 38.38 0.39512
## age 0.5195 1 39.58 0.47527
## air.temp 0.8094 1 1425.45 0.36845
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lmer.v1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: v1 ~ species + Sex + age + air.temp + (1 | Ring)
## Data: nat.var
##
## REML criterion at convergence: 1764.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.7165 -0.5902 0.0500 0.6787 3.4932
##
## Random effects:
## Groups Name Variance Std.Dev.
## Ring (Intercept) 0.05143 0.2268
## Residual 0.18525 0.4304
## Number of obs: 1434, groups: Ring, 49
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.808e+00 5.135e-01 3.940e+01 5.467 2.76e-06 ***
## speciesgreti 2.844e-01 9.757e-02 4.739e+01 2.914 0.00543 **
## speciesmarti 2.374e-01 2.101e-01 3.235e+01 1.130 0.26669
## SexF 7.792e-02 1.775e-01 3.166e+01 0.439 0.66366
## SexM 1.724e-01 1.762e-01 3.167e+01 0.978 0.33531
## age -6.299e-02 8.719e-02 3.862e+01 -0.722 0.47436
## air.temp -7.243e-03 8.038e-03 1.425e+03 -0.901 0.36768
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) spcsgr spcsmr SexF SexM age
## speciesgret -0.219
## speciesmart -0.444 0.049
## SexF -0.418 -0.057 0.763
## SexM -0.299 -0.149 0.742 0.895
## age -0.931 0.239 0.204 0.125 -0.004
## air.temp -0.161 -0.023 -0.008 0.003 -0.007 -0.008
lsmeans(lmer.v1, pairwise~species)
## $lsmeans
## species lsmean SE df lower.CL upper.CL
## bluti 2.47 0.0634 32.8 2.34 2.59
## greti 2.75 0.1061 44.1 2.54 2.96
## marti 2.70 0.1746 34.2 2.35 3.06
##
## Results are averaged over the levels of: Sex, age
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## bluti - greti -0.2844 0.0978 48.6 -2.906 0.0149
## bluti - marti -0.2374 0.2101 33.2 -1.130 0.5028
## greti - marti 0.0469 0.2273 35.2 0.207 0.9768
##
## Results are averaged over the levels of: Sex, age
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates
Anova(lmer.v2, type = 3, test.statistic = "F") #nothing significant
## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
##
## Response: v2
## F Df Df.res Pr(>F)
## (Intercept) 45.5217 1 41.03 3.741e-08 ***
## species 0.0849 2 40.16 0.9188
## Sex 0.0149 2 39.80 0.9852
## age 0.5283 1 40.47 0.4715
## air.temp 0.0280 1 1427.00 0.8671
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lmer.v2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: v2 ~ species + Sex + age + air.temp + (1 | Ring)
## Data: nat.var
##
## REML criterion at convergence: 1283
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.7896 -0.4106 0.0900 0.5407 3.6279
##
## Random effects:
## Groups Name Variance Std.Dev.
## Ring (Intercept) 0.04533 0.2129
## Residual 0.13154 0.3627
## Number of obs: 1434, groups: Ring, 49
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.202e+00 4.736e-01 3.859e+01 6.760 4.83e-08 ***
## speciesgreti -2.248e-02 8.931e-02 4.566e+01 -0.252 0.802
## speciesmarti -6.624e-02 1.955e-01 3.201e+01 -0.339 0.737
## SexF 1.907e-03 1.654e-01 3.152e+01 0.012 0.991
## SexM 1.414e-02 1.641e-01 3.152e+01 0.086 0.932
## age -5.862e-02 8.048e-02 3.806e+01 -0.728 0.471
## air.temp 1.139e-03 6.792e-03 1.427e+03 0.168 0.867
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) spcsgr spcsmr SexF SexM age
## speciesgret -0.221
## speciesmart -0.447 0.049
## SexF -0.420 -0.058 0.765
## SexM -0.307 -0.148 0.746 0.899
## age -0.932 0.240 0.203 0.123 0.000
## air.temp -0.148 -0.022 -0.007 0.003 -0.006 -0.008
Anova(lmer.avg.v, type = 3, test.statistic = "F")
## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
##
## Response: avg.v
## F Df Df.res Pr(>F)
## (Intercept) 42.9523 1 41.16 6.934e-08 ***
## species 1.1445 2 40.28 0.3285
## Sex 0.3973 2 40.09 0.6748
## age 0.7669 1 40.64 0.3863
## air.temp 0.1784 1 1426.98 0.6728
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lmer.avg.v)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: avg.v ~ species + Sex + age + air.temp + (1 | Ring)
## Data: nat.var
##
## REML criterion at convergence: 1137.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.5561 -0.4192 0.0955 0.6236 3.8317
##
## Random effects:
## Groups Name Variance Std.Dev.
## Ring (Intercept) 0.04272 0.2067
## Residual 0.11863 0.3444
## Number of obs: 1434, groups: Ring, 49
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.009e+00 4.582e-01 4.011e+01 6.566 7.48e-08 ***
## speciesgreti 1.253e-01 8.627e-02 4.729e+01 1.453 0.153
## speciesmarti 9.603e-02 1.895e-01 3.334e+01 0.507 0.616
## SexF 4.908e-02 1.603e-01 3.287e+01 0.306 0.761
## SexM 1.023e-01 1.591e-01 3.287e+01 0.643 0.525
## age -6.834e-02 7.787e-02 3.960e+01 -0.878 0.385
## air.temp -2.729e-03 6.454e-03 1.427e+03 -0.423 0.672
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) spcsgr spcsmr SexF SexM age
## speciesgret -0.221
## speciesmart -0.448 0.049
## SexF -0.421 -0.058 0.766
## SexM -0.308 -0.148 0.747 0.900
## age -0.933 0.240 0.203 0.122 0.000
## air.temp -0.145 -0.022 -0.007 0.003 -0.006 -0.008
Anova(lmer.acc, type = 3, test.statistic = "F")
## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
##
## Response: acc
## F Df Df.res Pr(>F)
## (Intercept) 3.1432 1 37.92 0.08428 .
## species 13.1702 2 37.23 4.728e-05 ***
## Sex 0.4682 2 33.21 0.63018
## age 0.6599 1 35.83 0.42195
## air.temp 0.8469 1 1374.24 0.35758
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lmer.acc)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: acc ~ species + Sex + age + air.temp + (1 | Ring)
## Data: nat.var
##
## REML criterion at convergence: 8283.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3928 -0.6556 -0.0021 0.6449 3.3918
##
## Random effects:
## Groups Name Variance Std.Dev.
## Ring (Intercept) 2.119 1.456
## Residual 18.195 4.266
## Number of obs: 1434, groups: Ring, 49
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 6.44546 3.62535 33.99118 1.778 0.0844 .
## speciesgreti -3.43432 0.71333 43.71998 -4.814 1.8e-05 ***
## speciesmarti -2.93246 1.42597 26.64274 -2.056 0.0497 *
## SexF -0.64160 1.19660 25.20720 -0.536 0.5965
## SexM -1.01659 1.18782 25.19098 -0.856 0.4002
## age -0.49848 0.61173 32.11682 -0.815 0.4211
## air.temp 0.07240 0.07845 1368.22941 0.923 0.3562
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) spcsgr spcsmr SexF SexM age
## speciesgret -0.201
## speciesmart -0.428 0.047
## SexF -0.406 -0.056 0.751
## SexM -0.266 -0.149 0.725 0.880
## age -0.923 0.223 0.208 0.132 -0.020
## air.temp -0.224 -0.024 -0.010 0.000 -0.011 -0.011
lsmeans(lmer.acc, pairwise~ species)
## $lsmeans
## species lsmean SE df lower.CL upper.CL
## bluti 3.946 0.429 28.7 3.07 4.82
## greti 0.511 0.763 43.5 -1.03 2.05
## marti 1.013 1.195 31.4 -1.42 3.45
##
## Results are averaged over the levels of: Sex, age
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## bluti - greti 3.434 0.716 48.7 4.795 <.0001
## bluti - marti 2.932 1.427 29.7 2.055 0.1167
## greti - marti -0.502 1.566 32.7 -0.320 0.9451
##
## Results are averaged over the levels of: Sex, age
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates
summary(model.avg(dredge(lmer(v1 ~ species + Sex + age + headwind + (1|Ring), data = nat.var, na.action = "na.fail")), subset = delta < 10))
## Warning in dredge(lmer(v1 ~ species + Sex + age + headwind + (1 | Ring), :
## comparing models fitted by REML
## Fixed term is "(Intercept)"
## Warning in dredge(lmer(v1 ~ species + Sex + age + headwind + (1 | Ring), :
## comparing models fitted by REML
## Fixed term is "(Intercept)"
##
## Call:
## model.avg(object = dredge(lmer(v1 ~ species + Sex + age + headwind +
## (1 | Ring), data = nat.var, na.action = "na.fail")), subset = delta <
## 10)
##
## Component model call:
## lmer(formula = v1 ~ <9 unique rhs>, data = nat.var, na.action =
## na.fail)
##
## Component models:
## df logLik AICc delta weight
## 4 5 -875.65 1761.34 0.00 0.55
## (Null) 3 -878.91 1763.84 2.50 0.16
## 24 6 -875.91 1763.87 2.53 0.16
## 14 6 -877.12 1766.30 4.96 0.05
## 2 4 -879.64 1767.31 5.96 0.03
## 1 4 -879.73 1767.50 6.15 0.03
## 34 7 -877.36 1768.80 7.46 0.01
## 124 7 -877.41 1768.89 7.55 0.01
## 12 5 -880.45 1770.95 9.60 0.00
##
## Term codes:
## age headwind Sex species
## 1 2 3 4
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 2.531459 0.188803 0.188894 13.401 <2e-16 ***
## speciesgreti 0.252291 0.155358 0.155394 1.624 0.104
## speciesmarti 0.113614 0.128321 0.128404 0.885 0.376
## headwind -0.006053 0.013160 0.013162 0.460 0.646
## age -0.005385 0.031687 0.031703 0.170 0.865
## SexF 0.001259 0.022706 0.022721 0.055 0.956
## SexM 0.002271 0.027944 0.027956 0.081 0.935
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 2.53146 0.18880 0.18889 13.401 < 2e-16 ***
## speciesgreti 0.32226 0.09101 0.09109 3.538 0.000403 ***
## speciesmarti 0.14512 0.12830 0.12840 1.130 0.258395
## headwind -0.02998 0.01185 0.01186 2.528 0.011466 *
## age -0.06034 0.08908 0.08914 0.677 0.498505
## SexF 0.09467 0.17298 0.17313 0.547 0.584480
## SexM 0.17075 0.17305 0.17320 0.986 0.324199
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model.avg(dredge(lmer(v2 ~ species + Sex + age + headwind + (1|Ring), data = nat.var, na.action = "na.fail")), subset = delta < 10))
## Warning in dredge(lmer(v2 ~ species + Sex + age + headwind + (1 | Ring), :
## comparing models fitted by REML
## Fixed term is "(Intercept)"
## Warning in dredge(lmer(v2 ~ species + Sex + age + headwind + (1 | Ring), :
## comparing models fitted by REML
## Fixed term is "(Intercept)"
##
## Call:
## model.avg(object = dredge(lmer(v2 ~ species + Sex + age + headwind +
## (1 | Ring), data = nat.var, na.action = "na.fail")), subset = delta <
## 10)
##
## Component model call:
## lmer(formula = v2 ~ <3 unique rhs>, data = nat.var, na.action =
## na.fail)
##
## Component models:
## df logLik AICc delta weight
## 2 4 -623.79 1255.61 0.00 0.92
## 12 5 -625.32 1260.68 5.07 0.07
## 23 6 -626.57 1265.20 9.59 0.01
##
## Term codes:
## age headwind species
## 1 2 3
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 2.775e+00 1.258e-01 1.258e-01 22.048 < 2e-16 ***
## headwind -4.589e-02 9.912e-03 9.920e-03 4.626 3.7e-06 ***
## age -3.349e-03 2.222e-02 2.224e-02 0.151 0.880
## speciesgreti 2.414e-05 7.129e-03 7.135e-03 0.003 0.997
## speciesmarti -3.536e-04 1.075e-02 1.075e-02 0.033 0.974
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 2.774741 0.125768 0.125848 22.048 < 2e-16 ***
## headwind -0.045892 0.009912 0.009920 4.626 3.7e-06 ***
## age -0.046010 0.069432 0.069491 0.662 0.508
## speciesgreti 0.003173 0.081676 0.081746 0.039 0.969
## speciesmarti -0.046479 0.114164 0.114261 0.407 0.684
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model.avg(dredge(lmer(avg.v ~ species + Sex + age + headwind + (1|Ring), data = nat.var, na.action = "na.fail")), subset = delta < 10))
## Warning in dredge(lmer(avg.v ~ species + Sex + age + headwind + (1 | Ring), :
## comparing models fitted by REML
## Fixed term is "(Intercept)"
## Warning in dredge(lmer(avg.v ~ species + Sex + age + headwind + (1 | Ring), :
## comparing models fitted by REML
## Fixed term is "(Intercept)"
##
## Call:
## model.avg(object = dredge(lmer(avg.v ~ species + Sex + age +
## headwind + (1 | Ring), data = nat.var, na.action = "na.fail")),
## subset = delta < 10)
##
## Component model call:
## lmer(formula = avg.v ~ <5 unique rhs>, data = nat.var, na.action =
## na.fail)
##
## Component models:
## df logLik AICc delta weight
## 2 4 -554.97 1117.97 0.00 0.82
## 12 5 -556.04 1122.13 4.16 0.10
## 24 6 -555.73 1123.52 5.54 0.05
## (Null) 3 -559.79 1125.59 7.62 0.02
## 23 6 -557.93 1127.93 9.96 0.01
##
## Term codes:
## age headwind Sex species
## 1 2 3 4
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 2.639e+00 1.811e-01 1.811e-01 14.569 < 2e-16 ***
## headwind -3.856e-02 1.073e-02 1.074e-02 3.591 0.000329 ***
## age -8.302e-03 3.299e-02 3.301e-02 0.252 0.801402
## speciesgreti 8.398e-03 4.021e-02 4.022e-02 0.209 0.834589
## speciesmarti 2.826e-03 2.765e-02 2.767e-02 0.102 0.918639
## SexF 5.814e-05 7.548e-03 7.554e-03 0.008 0.993859
## SexM 3.090e-04 8.433e-03 8.438e-03 0.037 0.970789
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 2.639152 0.181071 0.181144 14.569 < 2e-16 ***
## headwind -0.039275 0.009445 0.009453 4.155 3.26e-05 ***
## age -0.080695 0.068838 0.068897 1.171 0.241
## speciesgreti 0.163409 0.078297 0.078363 2.085 0.037 *
## speciesmarti 0.054993 0.109572 0.109665 0.501 0.616
## SexF 0.010274 0.099809 0.099893 0.103 0.918
## SexM 0.054606 0.097990 0.098073 0.557 0.578
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model.avg(dredge(lmer(acc ~ species + Sex + age + headwind + (1|Ring), data = nat.var, na.action = "na.fail")), subset = delta < 10))
## Warning in dredge(lmer(acc ~ species + Sex + age + headwind + (1 | Ring), :
## comparing models fitted by REML
## Fixed term is "(Intercept)"
## Warning in dredge(lmer(acc ~ species + Sex + age + headwind + (1 | Ring), :
## comparing models fitted by REML
## Fixed term is "(Intercept)"
##
## Call:
## model.avg(object = dredge(lmer(acc ~ species + Sex + age + headwind +
## (1 | Ring), data = nat.var, na.action = "na.fail")), subset = delta <
## 10)
##
## Component model call:
## lmer(formula = acc ~ <8 unique rhs>, data = nat.var, na.action =
## na.fail)
##
## Component models:
## df logLik AICc delta weight
## 34 7 -4141.33 8296.74 0.00 0.18
## 4 5 -4143.35 8296.75 0.01 0.18
## 14 6 -4142.43 8296.92 0.18 0.17
## 134 8 -4140.58 8297.26 0.52 0.14
## 234 8 -4140.97 8298.04 1.30 0.09
## 24 6 -4143.08 8298.21 1.47 0.09
## 124 7 -4142.17 8298.42 1.68 0.08
## 1234 9 -4140.25 8298.62 1.88 0.07
##
## Term codes:
## age headwind Sex species
## 1 2 3 4
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 5.28324 2.77238 2.77400 1.905 0.0568 .
## SexF -0.27265 0.88554 0.88621 0.308 0.7583
## SexM -0.49849 0.98070 0.98130 0.508 0.6115
## speciesgreti -3.40807 0.70124 0.70182 4.856 1.2e-06 ***
## speciesmarti -2.41490 1.23447 1.23543 1.955 0.0506 .
## age -0.25223 0.48361 0.48389 0.521 0.6022
## headwind -0.06771 0.11703 0.11706 0.578 0.5630
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 5.2832 2.7724 2.7740 1.905 0.0568 .
## SexF -0.5594 1.2036 1.2046 0.464 0.6424
## SexM -1.0228 1.1988 1.1998 0.852 0.3940
## speciesgreti -3.4081 0.7012 0.7018 4.856 1.2e-06 ***
## speciesmarti -2.4149 1.2345 1.2354 1.955 0.0506 .
## age -0.5537 0.5886 0.5891 0.940 0.3473
## headwind -0.2044 0.1158 0.1159 1.764 0.0778 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model.avg(dredge(lmer(v1 ~ species + Sex + age + air.temp + (1|Ring), data = nat.var, na.action = "na.fail")), subset = delta < 10))
## Warning in dredge(lmer(v1 ~ species + Sex + age + air.temp + (1 | Ring), :
## comparing models fitted by REML
## Fixed term is "(Intercept)"
## Warning in dredge(lmer(v1 ~ species + Sex + age + air.temp + (1 | Ring), :
## comparing models fitted by REML
## Fixed term is "(Intercept)"
##
## Call:
## model.avg(object = dredge(lmer(v1 ~ species + Sex + age + air.temp +
## (1 | Ring), data = nat.var, na.action = "na.fail")), subset = delta <
## 10)
##
## Component model call:
## lmer(formula = v1 ~ <6 unique rhs>, data = nat.var, na.action =
## na.fail)
##
## Component models:
## df logLik AICc delta weight
## 4 5 -875.65 1761.34 0.00 0.69
## (Null) 3 -878.91 1763.84 2.50 0.20
## 14 6 -877.12 1766.30 4.96 0.06
## 1 4 -879.73 1767.50 6.15 0.03
## 34 7 -877.36 1768.80 7.46 0.02
## 24 6 -879.18 1770.42 9.07 0.01
##
## Term codes:
## age air.temp Sex species
## 1 2 3 4
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 2.553e+00 1.861e-01 1.862e-01 13.711 <2e-16 ***
## speciesgreti 2.471e-01 1.569e-01 1.570e-01 1.574 0.115
## speciesmarti 1.117e-01 1.292e-01 1.293e-01 0.864 0.388
## age -5.519e-03 3.208e-02 3.210e-02 0.172 0.863
## SexF 1.566e-03 2.531e-02 2.533e-02 0.062 0.951
## SexM 2.824e-03 3.114e-02 3.115e-02 0.091 0.928
## air.temp -5.167e-05 9.139e-04 9.143e-04 0.057 0.955
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 2.552511 0.186079 0.186172 13.711 < 2e-16 ***
## speciesgreti 0.320657 0.091480 0.091557 3.502 0.000461 ***
## speciesmarti 0.145018 0.129793 0.129901 1.116 0.264260
## age -0.061645 0.089636 0.089703 0.687 0.491948
## SexF 0.094674 0.172979 0.173125 0.547 0.584480
## SexM 0.170751 0.173052 0.173199 0.986 0.324199
## air.temp -0.007001 0.008032 0.008039 0.871 0.383801
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model.avg(dredge(lmer(v2 ~ species + Sex + age + air.temp + (1|Ring), data = nat.var, na.action = "na.fail")), subset = delta < 10)) #null AICc higher by 11 delta AIC
## Warning in dredge(lmer(v2 ~ species + Sex + age + air.temp + (1 | Ring), :
## comparing models fitted by REML
## Fixed term is "(Intercept)"
## Warning in dredge(lmer(v2 ~ species + Sex + age + air.temp + (1 | Ring), :
## comparing models fitted by REML
## Fixed term is "(Intercept)"
##
## Call:
## model.avg(object = dredge(lmer(v2 ~ species + Sex + age + air.temp +
## (1 | Ring), data = nat.var, na.action = "na.fail")), subset = delta <
## 10)
##
## Component model call:
## lmer(formula = v2 ~ <3 unique rhs>, data = nat.var, na.action =
## na.fail)
##
## Component models:
## df logLik AICc delta weight
## (Null) 3 -630.73 1267.48 0.00 0.92
## 1 4 -632.26 1272.55 5.07 0.07
## 2 5 -633.45 1276.95 9.47 0.01
##
## Term codes:
## age species
## 1 2
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 2.918e+00 1.230e-01 1.230e-01 23.715 <2e-16 ***
## age -3.223e-03 2.239e-02 2.240e-02 0.144 0.886
## speciesgreti -5.486e-05 7.511e-03 7.517e-03 0.007 0.994
## speciesmarti -4.180e-04 1.147e-02 1.148e-02 0.036 0.971
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 2.917751 0.122954 0.123032 23.715 <2e-16 ***
## age -0.044292 0.071203 0.071263 0.622 0.534
## speciesgreti -0.006806 0.083375 0.083446 0.082 0.935
## speciesmarti -0.051853 0.116846 0.116945 0.443 0.657
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model.avg(dredge(lmer(avg.v ~ species + Sex + age + air.temp + (1|Ring), data = nat.var, na.action = "na.fail")), subset = delta < 10)) #null AICc higher by 7.62 delta AIC
## Warning in dredge(lmer(avg.v ~ species + Sex + age + air.temp + (1 | Ring), :
## comparing models fitted by REML
## Fixed term is "(Intercept)"
## Warning in dredge(lmer(avg.v ~ species + Sex + age + air.temp + (1 | Ring), :
## comparing models fitted by REML
## Fixed term is "(Intercept)"
##
## Call:
## model.avg(object = dredge(lmer(avg.v ~ species + Sex + age +
## air.temp + (1 | Ring), data = nat.var, na.action = "na.fail")),
## subset = delta < 10)
##
## Component model call:
## lmer(formula = avg.v ~ <4 unique rhs>, data = nat.var, na.action =
## na.fail)
##
## Component models:
## df logLik AICc delta weight
## (Null) 3 -559.79 1125.59 0.00 0.85
## 1 4 -560.88 1129.79 4.20 0.10
## 3 5 -560.83 1131.69 6.10 0.04
## 2 5 -562.63 1135.30 9.71 0.01
##
## Term codes:
## age Sex species
## 1 2 3
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 2.760e+00 1.809e-01 1.810e-01 15.256 <2e-16 ***
## age -8.309e-03 3.346e-02 3.347e-02 0.248 0.804
## speciesgreti 6.210e-03 3.440e-02 3.440e-02 0.181 0.857
## speciesmarti 2.033e-03 2.486e-02 2.487e-02 0.082 0.935
## SexF 4.366e-05 8.383e-03 8.390e-03 0.005 0.996
## SexM 3.882e-04 9.488e-03 9.494e-03 0.041 0.967
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 2.760498 0.180876 0.180950 15.256 <2e-16 ***
## age -0.079872 0.071031 0.071091 1.124 0.2612
## speciesgreti 0.154424 0.080822 0.080890 1.909 0.0563 .
## speciesmarti 0.050550 0.113628 0.113724 0.444 0.6567
## SexF 0.006592 0.102790 0.102877 0.064 0.9489
## SexM 0.058606 0.100899 0.100985 0.580 0.5617
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model.avg(dredge(lmer(acc ~ species + Sex + age + air.temp + (1|Ring), data = nat.var, na.action = "na.fail")), subset = delta < 10))
## Warning in dredge(lmer(acc ~ species + Sex + age + air.temp + (1 | Ring), :
## comparing models fitted by REML
## Fixed term is "(Intercept)"
## Warning in dredge(lmer(acc ~ species + Sex + age + air.temp + (1 | Ring), :
## comparing models fitted by REML
## Fixed term is "(Intercept)"
##
## Call:
## model.avg(object = dredge(lmer(acc ~ species + Sex + age + air.temp +
## (1 | Ring), data = nat.var, na.action = "na.fail")), subset = delta <
## 10)
##
## Component model call:
## lmer(formula = acc ~ <8 unique rhs>, data = nat.var, na.action =
## na.fail)
##
## Component models:
## df logLik AICc delta weight
## 34 7 -4141.33 8296.74 0.00 0.25
## 4 5 -4143.35 8296.75 0.01 0.24
## 14 6 -4142.43 8296.92 0.18 0.22
## 134 8 -4140.58 8297.26 0.52 0.19
## 234 8 -4142.53 8301.17 4.43 0.03
## 24 6 -4144.59 8301.24 4.50 0.03
## 124 7 -4143.66 8301.40 4.66 0.02
## 1234 9 -4141.78 8301.69 4.95 0.02
##
## Term codes:
## age air.temp Sex species
## 1 2 3 4
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 5.445073 2.765768 2.767390 1.968 0.0491 *
## SexF -0.274225 0.877139 0.877804 0.312 0.7547
## SexM -0.489946 0.970242 0.970840 0.505 0.6138
## speciesgreti -3.430691 0.697086 0.697663 4.917 9e-07 ***
## speciesmarti -2.421435 1.225460 1.226405 1.974 0.0483 *
## age -0.256582 0.484228 0.484500 0.530 0.5964
## air.temp 0.006876 0.032196 0.032211 0.213 0.8310
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 5.44507 2.76577 2.76739 1.968 0.0491 *
## SexF -0.56945 1.19564 1.19665 0.476 0.6342
## SexM -1.01741 1.19087 1.19188 0.854 0.3933
## speciesgreti -3.43069 0.69709 0.69766 4.917 9e-07 ***
## speciesmarti -2.42144 1.22546 1.22641 1.974 0.0483 *
## age -0.56042 0.58469 0.58518 0.958 0.3382
## air.temp 0.07084 0.07841 0.07848 0.903 0.3667
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lmer.v1.fin <- lmer(v1 ~ species + (1|Ring), data = nat.var)
lmer.v2.fin <- lmer(v2 ~ headwind + (1|Ring), data = nat.var)
lmer.avg.v.fin <- lmer(avg.v ~ headwind + (1|Ring), data = nat.var)
lmer.acc.fin <- lmer(acc ~ species + Sex + (1|Ring), data = nat.var)
#without outliers
lmer.v1.fin.out <- lmer(v1 ~ species + (1|Ring), data = nat.var.out)
lmer.v2.fin.out <- lmer(v2 ~ headwind + (1|Ring), data = nat.var.out)
lmer.avg.v.fin.out <- lmer(avg.v ~ headwind + (1|Ring), data = nat.var.out)
lmer.acc.fin.out <- lmer(acc ~ species + Sex + (1|Ring), data = nat.var.out)
Acceleration is skipped as the best model was the full model, acceleration only had one outlier, and the full model assumptions were well validated.
# Computes the default residuals (raw conditional)
raw_cond.v1.fin <- compute_redres(lmer.v1.fin)
raw_cond.v1.fin.out <- compute_redres(lmer.v1.fin.out)
raw_cond.v2.fin <- compute_redres(lmer.v2.fin)
raw_cond.v2.fin.out <- compute_redres(lmer.v2.fin.out)
raw_cond.avg.v.fin <- compute_redres(lmer.avg.v.fin)
raw_cond.avg.v.fin.out <- compute_redres(lmer.avg.v.fin.out)
#however checking the one without the outliers
raw_cond.acc.fin <- compute_redres(lmer.acc.fin)
raw_cond.acc.fin.out <- compute_redres(lmer.acc.fin.out)
shapiro.test(raw_cond.v1.fin)
##
## Shapiro-Wilk normality test
##
## data: raw_cond.v1.fin
## W = 0.96949, p-value < 2.2e-16
shapiro.test(raw_cond.v1.fin.out)
##
## Shapiro-Wilk normality test
##
## data: raw_cond.v1.fin.out
## W = 0.98207, p-value = 2.746e-12
shapiro.test(raw_cond.v2.fin)
##
## Shapiro-Wilk normality test
##
## data: raw_cond.v2.fin
## W = 0.9004, p-value < 2.2e-16
shapiro.test(raw_cond.v2.fin.out)
##
## Shapiro-Wilk normality test
##
## data: raw_cond.v2.fin.out
## W = 0.96104, p-value < 2.2e-16
shapiro.test(raw_cond.avg.v.fin)
##
## Shapiro-Wilk normality test
##
## data: raw_cond.avg.v.fin
## W = 0.90085, p-value < 2.2e-16
shapiro.test(raw_cond.avg.v.fin.out)
##
## Shapiro-Wilk normality test
##
## data: raw_cond.avg.v.fin.out
## W = 0.94798, p-value < 2.2e-16
shapiro.test(raw_cond.acc.fin)
##
## Shapiro-Wilk normality test
##
## data: raw_cond.acc.fin
## W = 0.9936, p-value = 7.438e-06
shapiro.test(raw_cond.acc.fin.out)
##
## Shapiro-Wilk normality test
##
## data: raw_cond.acc.fin.out
## W = 0.99372, p-value = 1.043e-05
shapiro.test(raw_cond.acc.fin)
##
## Shapiro-Wilk normality test
##
## data: raw_cond.acc.fin
## W = 0.9936, p-value = 7.438e-06
shapiro.test(raw_cond.acc.fin.out)
##
## Shapiro-Wilk normality test
##
## data: raw_cond.acc.fin.out
## W = 0.99372, p-value = 1.043e-05
hist(raw_cond.v1.fin) #normality in histogram
hist(raw_cond.v1.fin.out)
plot_redres(lmer.v1.fin) #heteroskedasticity
plot_redres(lmer.v1.fin.out)
plot_resqq(lmer.v1.fin) #normality of residuals
plot_resqq(lmer.v1.fin.out)
plot_ranef(lmer.v1.fin) #normality of random effects residuals
plot_ranef(lmer.v1.fin.out)
hist(raw_cond.v2.fin) #normality in histogram
hist(raw_cond.v2.fin.out)
plot_redres(lmer.v2.fin) #heteroskedasticity
plot_redres(lmer.v2.fin.out)
plot_resqq(lmer.v2.fin) #normality of residuals
plot_resqq(lmer.v2.fin.out)
plot_ranef(lmer.v2.fin) #normality of random effects residuals
plot_ranef(lmer.v2.fin.out)
hist(raw_cond.avg.v.fin) #normality in histogram
hist(raw_cond.avg.v.fin.out)
plot_redres(lmer.avg.v.fin) #heteroskedasticity
plot_redres(lmer.avg.v.fin.out)
plot_resqq(lmer.avg.v.fin) #normality of residuals
plot_resqq(lmer.avg.v.fin.out)
plot_ranef(lmer.avg.v.fin) #normality of random effects residuals
plot_ranef(lmer.avg.v.fin.out)
#contrasts = list(Ring = contr.sum, cgt_day = contr.sum))
Anova(lmer.v1.fin, type = 3, test.statistic = "F") #species
## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
##
## Response: v1
## F Df Df.res Pr(>F)
## (Intercept) 3494.6627 1 39.858 < 2.2e-16 ***
## species 6.3849 2 45.683 0.003588 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(lmer.v1.fin.out, type = 3, test.statistic = "F")
## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
##
## Response: v1
## F Df Df.res Pr(>F)
## (Intercept) 3720.5816 1 39.989 < 2.2e-16 ***
## species 6.1535 2 45.552 0.004309 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lmer.v1.fin)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: v1 ~ species + (1 | Ring)
## Data: nat.var
##
## REML criterion at convergence: 1751.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.7220 -0.5874 0.0470 0.6827 3.5530
##
## Random effects:
## Groups Name Variance Std.Dev.
## Ring (Intercept) 0.04864 0.2205
## Residual 0.18537 0.4305
## Number of obs: 1434, groups: Ring, 49
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.50680 0.04235 38.68447 59.190 < 2e-16 ***
## speciesgreti 0.32179 0.09112 52.30977 3.531 0.000872 ***
## speciesmarti 0.14344 0.12596 38.05612 1.139 0.261922
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) spcsgr
## speciesgret -0.465
## speciesmart -0.336 0.156
tab_model(lmer.v1.fin)
| v1 | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 2.51 | 2.42 – 2.59 | <0.001 |
| species [greti] | 0.32 | 0.14 – 0.50 | <0.001 |
| species [marti] | 0.14 | -0.10 – 0.39 | 0.255 |
| Random Effects | |||
| σ2 | 0.19 | ||
| τ00 Ring | 0.05 | ||
| ICC | 0.21 | ||
| N Ring | 49 | ||
| Observations | 1434 | ||
| Marginal R2 / Conditional R2 | 0.042 / 0.241 | ||
lsmeans(lmer.v1.fin, pairwise~species)
## $lsmeans
## species lsmean SE df lower.CL upper.CL
## bluti 2.51 0.0424 39.9 2.42 2.59
## greti 2.83 0.0809 59.0 2.67 2.99
## marti 2.65 0.1187 39.1 2.41 2.89
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## bluti - greti -0.322 0.0913 53.9 -3.524 0.0025
## bluti - marti -0.143 0.1260 39.2 -1.138 0.4967
## greti - marti 0.178 0.1436 44.2 1.242 0.4352
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates
lsmeans.v1 <- lsmeans(lmer.v1.fin, pairwise~species)
Anova(lmer.v2.fin, type = 3, test.statistic = "F") #species and headwind
## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
##
## Response: v2
## F Df Df.res Pr(>F)
## (Intercept) 3857.990 1 160.93 < 2.2e-16 ***
## headwind 21.387 1 1431.16 4.091e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(lmer.v2.fin.out, type = 3, test.statistic = "F")
## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
##
## Response: v2
## F Df Df.res Pr(>F)
## (Intercept) 4472.346 1 142.06 < 2.2e-16 ***
## headwind 21.269 1 1413.59 4.351e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lmer.v2.fin)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: v2 ~ headwind + (1 | Ring)
## Data: nat.var
##
## REML criterion at convergence: 1247.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.8745 -0.4077 0.0902 0.5318 3.5448
##
## Random effects:
## Groups Name Variance Std.Dev.
## Ring (Intercept) 0.03895 0.1973
## Residual 0.12968 0.3601
## Number of obs: 1434, groups: Ring, 49
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.757e+00 4.433e-02 1.571e+02 62.19 < 2e-16 ***
## headwind -4.589e-02 9.912e-03 1.431e+03 -4.63 3.99e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## headwind 0.700
Anova(lmer.avg.v.fin, type = 3, test.statistic = "F") #species and headwind
## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
##
## Response: avg.v
## F Df Df.res Pr(>F)
## (Intercept) 3587.828 1 149.85 < 2.2e-16 ***
## headwind 17.223 1 1429.98 3.52e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(lmer.avg.v.fin.out, type = 3, test.statistic = "F")
## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
##
## Response: avg.v
## F Df Df.res Pr(>F)
## (Intercept) 4221.690 1 135.54 < 2.2e-16 ***
## headwind 13.393 1 1412.39 0.0002619 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lmer.avg.v.fin)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: avg.v ~ headwind + (1 | Ring)
## Data: nat.var
##
## REML criterion at convergence: 1109.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.6928 -0.4278 0.1154 0.6222 3.7496
##
## Random effects:
## Groups Name Variance Std.Dev.
## Ring (Intercept) 0.03948 0.1987
## Residual 0.11744 0.3427
## Number of obs: 1434, groups: Ring, 49
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.595e+00 4.328e-02 1.428e+02 59.969 < 2e-16 ***
## headwind -3.924e-02 9.444e-03 1.430e+03 -4.155 3.45e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## headwind 0.683
Anova(lmer.acc.fin, type = 3, test.statistic = "F") #species and Sex
## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
##
## Response: acc
## F Df Df.res Pr(>F)
## (Intercept) 16.0819 1 28.151 0.0004057 ***
## species 12.7559 2 39.110 5.436e-05 ***
## Sex 0.6363 2 33.651 0.5354960
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(lmer.acc.fin.out, type = 3, test.statistic = "F")
## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
##
## Response: acc
## F Df Df.res Pr(>F)
## (Intercept) 15.522 1 28.253 0.0004878 ***
## species 13.214 2 39.044 4.141e-05 ***
## Sex 0.493 2 33.764 0.6150967
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lmer.acc.fin)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: acc ~ species + Sex + (1 | Ring)
## Data: nat.var
##
## REML criterion at convergence: 8282.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3766 -0.6433 0.0123 0.6439 3.4045
##
## Random effects:
## Groups Name Variance Std.Dev.
## Ring (Intercept) 2.162 1.470
## Residual 18.181 4.264
## Number of obs: 1434, groups: Ring, 49
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.4944 1.1206 27.0775 4.011 0.000429 ***
## speciesgreti -3.2851 0.6995 50.4970 -4.696 2.07e-05 ***
## speciesmarti -2.6814 1.4065 29.5907 -1.906 0.066347 .
## SexF -0.5125 1.1964 27.9608 -0.428 0.671682
## SexM -1.0269 1.1977 28.0081 -0.857 0.398506
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) spcsgr spcsmr SexF
## speciesgret 0.000
## speciesmart -0.797 0.000
## SexF -0.937 -0.089 0.746
## SexM -0.936 -0.149 0.745 0.891
tab_model(lmer.acc.fin)
| acc | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 4.49 | 2.30 – 6.69 | <0.001 |
| species [greti] | -3.29 | -4.66 – -1.91 | <0.001 |
| species [marti] | -2.68 | -5.44 – 0.08 | 0.057 |
| Sex [F] | -0.51 | -2.86 – 1.83 | 0.668 |
| Sex [M] | -1.03 | -3.38 – 1.32 | 0.391 |
| Random Effects | |||
| σ2 | 18.18 | ||
| τ00 Ring | 2.16 | ||
| ICC | 0.11 | ||
| N Ring | 49 | ||
| Observations | 1434 | ||
| Marginal R2 / Conditional R2 | 0.062 / 0.161 | ||
lsmeans(lmer.acc.fin, pairwise~species)
## $lsmeans
## species lsmean SE df lower.CL upper.CL
## bluti 3.981 0.429 29.6 3.105 4.86
## greti 0.696 0.737 46.9 -0.787 2.18
## marti 1.300 1.151 32.4 -1.044 3.64
##
## Results are averaged over the levels of: Sex
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## bluti - greti 3.285 0.702 52.5 4.679 0.0001
## bluti - marti 2.681 1.407 30.8 1.906 0.1541
## greti - marti -0.604 1.572 33.8 -0.384 0.9221
##
## Results are averaged over the levels of: Sex
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates
lsmeans.acc <- lsmeans(lmer.acc.fin, pairwise~species)
plot.v1.means <- data.frame(lsmeans.v1$lsmeans)
plot.acc.means <- data.frame(lsmeans.acc$lsmeans)
plot.v1.lsmeans <-
ggplot() +
geom_point(data = plot.v1.means, aes(x = species, y = lsmean, color = species), alpha = 1) +
geom_errorbar(data = plot.v1.means, aes(x = species, ymin = lower.CL, ymax = upper.CL, color = species, width = 0.5),alpha = 1) +
scale_color_manual(values = c("bluti" ="black", "greti" = "black", "marti" = "black")) +
scale_x_discrete(labels = c("Blue Tit","Great Tit","Marsh Tit")) +
ylab(expression(initial~velocity~(m/s))) +
ylim(2.25,3.25)+
theme(axis.title.x=element_blank()) + #getting rid of x-axis title
theme(legend.position = "none") #getting rid of legend
plot.acc.lsmeans <-
ggplot() +
geom_point(data = plot.acc.means, aes(x = species, y = lsmean, color = species), alpha = 1) +
geom_errorbar(data = plot.acc.means, aes(x = species, ymin = lower.CL, ymax = upper.CL, color = species, width = 0.5),alpha = 1) +
scale_color_manual(values = c("bluti" ="black", "greti" = "black", "marti" = "black")) +
scale_x_discrete(labels = c("Blue Tit","Great Tit","Marsh Tit")) +
ylab(expression(`in`-flight~acceleration~(m/s^2))) +
ylim(-2,6) +
theme(axis.title.x=element_blank()) + #getting rid of x-axis title
theme(legend.position = "none") #getting rid of legend
grid.arrange(plot.v1.lsmeans, plot.acc.lsmeans, nrow = 1)